All first-order linear ordinary differential equations are soluble, and simply enough that the method appears in all textbooks on differential equations. Textbooks commonly state that there is no corresponding general solution for second-order equations when in fact there is. It appeared for the first time in Tricomi’s Differential Equations and was reprised in his Integral Equations, but to the best of my knowledge appears in no other texts. Presumably that is because the method requires a basic understanding of integral equations.

Consider a general second-order linear ordinary differential equation with nonconstant coefficients:

$A\left(x\right)y^{″}+B\left(x\right)y^{\prime }+C\left(x\right)y=0$

The coefficient functions can be split additively

$\left[a_0\right(x)+a(x\left)\right]y^{″}+\left[b_0\right(x)+b(x\left)\right]y^{\prime }+\left[c_0\right(x)+c(x\left)\right]y=0$

in an infinite number of ways. The essence of Fubini’s method is to perform this split so that one set of coefficients satisfies an equation that can be solved explicitly,

$a_0\left(x\right)u^{″}+b_0\left(x\right)u^{\prime }+c_0\left(x\right)u=0\to u\left(x\right)=f_1\left(x\right),f_2\left(x\right)$

and then look for a solution to the original equation of the form

$y\left(x\right)=c_1\left(x\right)f_1\left(x\right)+c_2\left(x\right)f_2\left(x\right)$

along with the same constraint as is used in variation of parameters:

$c_1^{\prime }\left(x\right)f_1\left(x\right)+c_2^{\prime }\left(x\right)f_2\left(x\right)=0$

This form is now to be substituted in the original equation, where in the process of forming derivatives the imposed constraint removes half of the terms of the first derivative. The resulting equation is

$c_1^{\prime }{f}_1^{\prime }+c_2^{\prime }{f}_2^{\prime }+c_1\frac{a{f}_1^{″}+b{f}_1^{\prime }+c{f}_1}{a_0+a}+c_2\frac{a{f}_2^{″}+b{f}_2^{\prime }+c{f}_2}{a_0+a}=0$

where it is assumed that denominators remain nonzero. These two equations are easily solved for the derivatives of the varied parameters,

$c_1^{\prime }=(c_1{g}_1+c_2{g}_2)f_2{c}_2^{\prime }=-(c_1{g}_1+c_2{g}_2)f_1$

which for conciseness uses the abbreviations

$\begin{array}{c}{g}_1\left(x\right)=\frac{a{f}_1^{″}+b{f}_1^{\prime }+c{f}_1}{(a_0+a)W}{g}_2\left(x\right)=\frac{a{f}_2^{″}+b{f}_2^{\prime }+c{f}_2}{(a_0+a)W}\\ \\ W\left(x\right)=f_1{f}_2^{\prime }-f_1^{\prime }{f}_2\end{array}$

The Wronskian is nonzero by definition for the two linearly independent solutions of the intermediate differential equation. The varied parameters thus satisfy the system of Volterra integral equations

$\begin{array}{l}{c}_1\left(x\right)={\gamma }_1+\underset{0}{\overset{x}{\int }}dt\left[c_1\right(t\left)g_1\right(t)+c_2(t\left)g_2\right(t\left)\right]f_2\left(t\right)\\ c_2\left(x\right)={\gamma }_2-\underset{0}{\overset{x}{\int }}dt\left[c_1\right(t\left)g_1\right(t)+c_2(t\left)g_2\right(t\left)\right]f_1\left(t\right)\end{array}$

where the constants of integration ${\gamma }_i=c_i\left(0\right)$ are related to the initial conditions of the original equation by

$\begin{array}{c}{\gamma }_1{f}_1\left(0\right)+{\gamma }_2{f}_2\left(0\right)=y\left(0\right)\\ {\gamma }_1{f}_1^{\prime }\left(0\right)+{\gamma }_2{f}_2^{\prime }\left(0\right)=y^{\prime }\left(0\right)\\ \\ {\gamma }_1=\frac{f_2^{\prime }\left(0\right)y\left(0\right)-f_2\left(0\right)y^{\prime }\left(0\right)}{W\left(0\right)}{\gamma }_2=\frac{-f_1^{\prime }\left(0\right)y\left(0\right)+f_1\left(0\right)y^{\prime }\left(0\right)}{W\left(0\right)}\end{array}$

In order to decouple the integral equations, first set

$h\left(x\right)=c_1{g}_1+c_2{g}_2$

for the simpler equations

$c_1\left(x\right)={\gamma }_1+\underset{0}{\overset{x}{\int }}dth\left(t\right)f_2\left(t\right)c_2\left(x\right)={\gamma }_2-\underset{0}{\overset{x}{\int }}dth\left(t\right)f_1\left(t\right)$

Multiplying the first by g₁ and the second by g₂ then adding gives

$h\left(x\right)={\gamma }_1{g}_1\left(x\right)+{\gamma }_2{g}_2\left(x\right)+\underset{0}{\overset{x}{\int }}dt\left[g_1\right(x\left)f_2\right(t)-g_2(x\left)f_1\right(t\left)\right]h\left(t\right)$

When this equation is solved using successive approximations, the result will be linear in the constants of integration. That means one can assume

$h\left(x\right)={\gamma }_1{h}_1\left(x\right)+{\gamma }_2{h}_2\left(x\right)$

and arrive at the pair of decoupled Volterra integral equations:

$\begin{array}{c}{h}_i\left(x\right)=g_i\left(x\right)+\underset{0}{\overset{x}{\int }}dtK(x,t)h_i\left(t\right)i=1,2\\ \\ K(x,t)=g_1\left(x\right)f_2\left(t\right)-g_2\left(x\right)f_1\left(t\right)\end{array}$

To solve these equations by successive approximations, start each solution with the first term on the right-hand side and repeatedly substitute the result back into the equations:

$\begin{array}{l}{h}_{i,0}\left(x\right)=g_i\left(x\right)\\ h_{i,1}\left(x\right)=g_i\left(x\right)+\underset{0}{\overset{x}{\int }}dtK(x,t)g_i\left(t\right)\\ h_{i,2}\left(x\right)=g_i\left(x\right)+\underset{0}{\overset{x}{\int }}dtK(x,t)g_i\left(t\right)+\underset{0}{\overset{x}{\int }}dtK(x,t)\underset{0}{\overset{t}{\int }}duK(x,u)g_i\left(u\right)\end{array}$

Reversing the order of integration step by step, one can define the iterated kernels

$\begin{array}{l}{K}_1(x,t)=K(x,t)\\ K_n(x,t)=\underset{t}{\overset{x}{\int }}duK(x,u)K_{n-1}(u,t)\end{array}$

where the change in the lower limit of integration comes from the reversal. The solutions to the integral equations are then

$h_i\left(x\right)=g_i\left(x\right)+\sum _{n=1}^{\infty }\underset{0}{\overset{x}{\int }}dt{K}_n(x,t)g_i\left(t\right)i=1,2$

When the iterated kernels can be recognized as converging to a known function this solution has only a single integral. With these solutions substituted in the equations for the varied parameters, the general solution of the original differential equation, in terms of the intermediate constants of integration, is

$\begin{array}{c}y\left(x\right)={\gamma }_1{y}_1\left(x\right)+{\gamma }_2{y}_2\left(x\right)\\ \\ y_i\left(x\right)=f_i\left(x\right)+\underset{0}{\overset{x}{\int }}dt\left[f_1\right(x\left)f_2\right(t)-f_2(x\left)f_1\right(t\left)\right]h_i\left(t\right)\end{array}$

This result can easily be rewritten in terms of the initial conditions of the original differential equation, but this form is patently simpler.

Since an infinite sum appears in the integral equations determining h_i it might not appear that this method is an improvement over a power series solution to the original equation. This method at least provides a straightforward algorithm to the solutions, while a recurrence relation derived from a power series might not be apparently soluble. It is also easier in general to establish if this form converges to a solution for given coefficient functions in the original equation.

As a first example of Fubini’s method, consider the differential equation for hyperbolic circular functions

$y^{″}-y=0$

which will be solved with respect to the auxiliary equation

$u^{″}=0\to f_1=1,f_2=x$

The splitting coefficients of the original equation are

$\begin{array}{l}{a}_0=1\\ b_0=0\\ c_0=0\end{array}\begin{array}{l}a=0\\ b=0\\ c=-1\end{array}$

and the intermediate functions are

$W\left(x\right)=1\begin{array}{l}{g}_1\left(x\right)=-1\\ g_2\left(x\right)=-x\end{array}$

The first few iterated kernels are

$K_1(x,t)=K(x,t)=x-t$

$K_2(x,t)=\underset{t}{\overset{x}{\int }}du(x-u)(u-t)=\frac{1}{3!}(x-t{)}^3$

$K_3(x,t)=\frac{1}{3!}\underset{t}{\overset{x}{\int }}du(x-u)(u-t{)}^3=\frac{1}{5!}(x-t{)}^5$

$K_4(x,t)=\frac{1}{5!}\underset{t}{\overset{x}{\int }}du(x-u)(u-t{)}^5=\frac{1}{7!}(x-t{)}^7$

and the pattern is clear. In this case the iterated kernels can be summed to the single function

$\sum _{n=1}^{\infty }{K}_n(x,t)=sinh(x-t)$

The solutions to the decoupled integral equations are

$\begin{array}{l}{h}_1\left(x\right)=-1-\underset{0}{\overset{x}{\int }}dtsinh(x-t)=-cosh\left(x\right)\\ h_2\left(x\right)=-x-\underset{0}{\overset{x}{\int }}dttsinh(x-t)=-sinh\left(x\right)\end{array}$

and the solutions to the original differential equation are

$\begin{array}{l}{y}_1\left(x\right)=1+\underset{0}{\overset{x}{\int }}dt(x-t)cosh\left(t\right)=cosh\left(x\right)\\ y_2\left(x\right)=x+\underset{0}{\overset{x}{\int }}dt(x-t)sinh\left(t\right)=sinh\left(x\right)\end{array}$

which are precisely the hyperbolic circular functions as expected. Nifty!

Repeat the process for the Airy differential equation

$y^{″}-xy=0$

using the same auxiliary equation and solutions. The splitting coefficients of the original equation are

$\begin{array}{l}{a}_0=1\\ b_0=0\\ c_0=0\end{array}\begin{array}{l}a=0\\ b=0\\ c=-x\end{array}$

and the intermediate functions are

$W\left(x\right)=1\begin{array}{l}{g}_1\left(x\right)=-x\\ g_2\left(x\right)=-x^2\end{array}$

The first few iterated kernels are

$K_1(x,t)=K(x,t)=x(x-t)$

$K_2(x,t)=x\underset{t}{\overset{x}{\int }}du(x-u)u(u-t)=\frac{x}{12}(x-t{)}^3(x+t)$

$\begin{array}{l}{K}_3(x,t)=\frac{x}{12}\underset{t}{\overset{x}{\int }}du(x-u)u(u-t{)}^3(u+t)\\ \phantom{K_3(x,t)}=\frac{x}{2520}(x-t{)}^5(5x^2+11xt+5t^2)\end{array}$

$\begin{array}{l}{K}_4(x,t)=\frac{x}{2520}\underset{t}{\overset{x}{\int }}du(x-u)u(u-t{)}^3(5u^2+11ut+5t^2)\\ \phantom{K_4(x,t)}=\frac{x}{90720}(x-t{)}^7(x+t\left)\right(x+2t\left)\right(2x+t\left)\right)\end{array}$

This time the pattern is far from clear, but one can still find a solution to this order of accuracy. The solutions to the decoupled integral equations are

$\begin{array}{l}{h}_1\left(x\right)=-x-\underset{0}{\overset{x}{\int }}dtt(K_1+K_2+K_3+K_4+\cdots )\\ \phantom{h_1\left(x\right)}=-x-\frac{x^4}{6}-\frac{x^7}{180}-\frac{x^{10}}{12960}-\frac{x^{13}}{1710720}-\cdots \\ \\ h_2\left(x\right)=-x^2-\underset{0}{\overset{x}{\int }}dt{t}^2(K_1+K_2+K_3+K_4+\cdots )\\ \phantom{h_2\left(x\right)}=-x^2-\frac{x^5}{12}-\frac{x^8}{504}-\frac{x^{11}}{45360}-\frac{x^{14}}{7076160}-\cdots \end{array}$

and the two solutions to the original differential equation are

$\begin{array}{l}{y}_1\left(x\right)=1-\underset{0}{\overset{x}{\int }}dt(x-t)h_1\left(t\right)\\ \phantom{y_1\left(x\right)}=1+\frac{x^3}{6}+\frac{x^6}{180}+\frac{x^9}{12960}+\frac{x^{12}}{1710720}+\cdots \\ \\ y_2\left(x\right)=x-\underset{0}{\overset{x}{\int }}dt(x-t)h_2\left(t\right)\\ \phantom{y_2\left(x\right)}=x+\frac{x^4}{12}+\frac{x^7}{504}+\frac{x^{10}}{45360}+\frac{x^{13}}{7076160}+\cdots \\ \phantom{y_1\left(x\right)}\end{array}$

which curiously have the same numeric denominators as the solutions to the decoupled integral equations. Using power series expansions of combinations of Airy functions, one can make the identifications

$\begin{array}{l}{y}_1\left(x\right)=\frac{1}{2}[\frac{Ai\left(x\right)}{Ai\left(0\right)}+\frac{Bi\left(x\right)}{Bi\left(0\right)}]\\ y_2\left(x\right)=\frac{\Gamma \left(\frac{1}{3}\right)}{2·3^{1/3}\Gamma \left(\frac{2}{3}\right)}[\frac{Bi\left(x\right)}{Bi\left(0\right)}-\frac{Ai\left(x\right)}{Ai\left(0\right)}]\end{array}$

The solutions to the decoupled integral equations are clearly related to the final solutions by

$h_1\left(x\right)=-x{y}_1\left(x\right)h_2\left(x\right)=-x{y}_2\left(x\right)$

The equal numeric denominators can be understood term by term with the integral

$\underset{0}{\overset{x}{\int }}dt{t}^n(x-t)=\frac{x^{n+2}}{(n+1)(n+2)}$

The relation of the iterated kernels to the terms in the two Airy series remains to be explored.

A Mathematica notebook containing evaluation of all integrals can be found here.

Since any second-order linear ordinary differential equation can be transformed to the normal form

$y^{″}+f\left(x\right)y=0$

one can now write a general solution to this equation, using the same auxiliary equation and similar intermediate functions as in the two examples. The iterated kernels are in general

$\begin{array}{l}{K}_1(x,t)=-f\left(x\right)(x-t)\\ K_n(x,t)=-f\left(x\right)\underset{t}{\overset{x}{\int }}du(x-u)K_{n-1}(u,t)\end{array}$

Once these are determined, evaluate the functions

$\begin{array}{l}{h}_1\left(x\right)=f\left(x\right)+\sum _{n=1}^{\infty }\underset{0}{\overset{x}{\int }}dt{K}_n(x,t)f\left(t\right)\\ h_2\left(x\right)=xf\left(x\right)+\sum _{n=1}^{\infty }\underset{0}{\overset{x}{\int }}dt{K}_n(x,t)tf\left(t\right)\end{array}$

Then the solutions to the differential equation are simply

$y_1\left(x\right)=1-\underset{0}{\overset{x}{\int }}dt(x-t)h_1\left(t\right)y_2\left(x\right)=x-\underset{0}{\overset{x}{\int }}dt(x-t)h_2\left(t\right)$

The main difficulty is of course the evaluation of the iterated kernels. This is again a completely straightforward algorithm to the general solution. As such it deserves to be much better known than currently.

Uploaded 2020.08.25 analyticphysics.com